module adv_1d_case_mod
  
  use constants_mod
  use namelist_mod
  use boundary_mod, only: set_periodic_boundary

  implicit none

  private
  public :: adv_1d_case_init
  public :: get_exact_solution

contains
    
  subroutine adv_1d_case_init(q)

    real(r8), intent(out) :: q(1-nhalo:nx+nhalo)
    integer :: i

    select case (trim(test_case))
    case ("sine_wave")
      do i = 1, nx
        q(i) = sin(pi * x(i))
      end do
    case ('cosine')
      do i = 1, nx
        if (x(i) >= 0.25_r8 .and. x(i) <= 0.35_r8) then
          q(i) = 1.0_r8 + cos(pi * (x(i) - 0.3_r8) / 0.05_r8)
        else
          q(i) = 0.0_r8
        end if 
      end do
    case ('gaussian')
      do i = 1, nx
        q(i) = 1.0_r8 / (1.0_r8 + exp(80.0_r8 * (abs(x(i)) - 0.15_r8)))
      end do
    case ('square_wave')
      do i = 1, nx
        if (abs(x(i)) <= 0.25_r8) then
        ! if (x(i) <= 30.0_r8) then
          q(i) = 1.0_r8
        else
          q(i) = 0.0_r8
        end if
      end do
    case default
      write(*,*) 'Unknown test case: ', trim(test_case)
      stop
    end select
    
    call set_periodic_boundary(q)

  end subroutine adv_1d_case_init

  subroutine get_exact_solution(elaplsed_time, rho_exact)

    real(r8), intent(in) :: elaplsed_time
    real(r8), intent(out) :: rho_exact(1:nx)
    integer i, n_periods
    real(r8) shifted_x1, shifted_x2, L, x1, x2
    logical in_wave
    real(r8) total_shift, wave_center, wave_left, wave_right, shifted_center, shifted_left, shifted_right

    L = x_max - x_min

    select case (trim(test_case))
    case ('sine_wave')
      do i = 1, nx
        shifted_x1 = x(i) - u * elaplsed_time
        shifted_x1 = mod(shifted_x1, L)
        rho_exact(i) = sin(pi * shifted_x1)
      end do 

    case ('square_wave')
      total_shift = u * elaplsed_time
      ! 方波初始中心在0，宽度为0.5 [-0.25, 0.25]
      wave_center = 0.0_r8
      wave_left = -0.25_r8
      wave_right = 0.25_r8
      ! 计算平流后的位置
      shifted_center = wave_center + total_shift
      shifted_left = wave_left + total_shift
      shifted_right = wave_right + total_shift
      
      do i = 1, nx
        ! 关键：使用周期性映射，确保位置在 [x_min, x_max] 内
        call apply_periodic_bc(x_min, x_max, shifted_left, shifted_right, x(i), in_wave)
        rho_exact(i) = merge(1.0_r8, 0.0_r8, in_wave)
      end do
    
    end select

  end subroutine get_exact_solution

  subroutine apply_periodic_bc(x_min, x_max, wave_left, wave_right, x, in_wave)

    real(r8), intent(in) :: x_min, x_max, wave_left, wave_right, x
    logical, intent(out) :: in_wave
    
    real(r8) :: L, left_mod, right_mod
    real(r8) :: wave_width
    
    L = x_max - x_min
    wave_width = wave_right - wave_left
    
    ! 将波的位置映射到基本区间 [x_min, x_max]
    left_mod = modulo(wave_left - x_min, L) + x_min
    right_mod = left_mod + wave_width  ! 保持波宽度不变
    
    ! 处理跨越边界的情况
    if (right_mod <= x_max) then
      ! 波完全在域内
      in_wave = (x >= left_mod .and. x <= right_mod)
    else
      ! 波跨越右边界
      in_wave = (x >= left_mod .and. x <= x_max) .or. &
                (x >= x_min .and. x <= right_mod - L)
    end if
    
end subroutine apply_periodic_bc

end module adv_1d_case_mod